# Title: script used in "Clientelism and Local Politics: Interactions between Councilors and Voters in the State of Minas Gerais"

# Date: 18/07/2019
#Libraries
library(ordinal)
library(tidyverse)
library(plyr)
library(lattice)

#reading tha data
dataset<-readRDS("Clientelism_data.rds")

#or

dataset <- read.csv2("Clientelism_data.csv")

#Table 1 Chi-squared test: table 1
#####
with(dataset,chisq.test(V56rec,V57Arec))
with(dataset,chisq.test(V56rec,V58Arec))
with(dataset,chisq.test(V56rec,V60Arec))
with(dataset,chisq.test(V56rec,V61Arec))
with(dataset,chisq.test(V57Arec,V58Arec))
with(dataset,chisq.test(V57Arec,V60Arec))
with(dataset,chisq.test(V57Arec,V61Arec))
with(dataset,chisq.test(V58Arec,V60Arec))
with(dataset,chisq.test(V58Arec,V61Arec))
with(dataset,chisq.test(V60Arec,V61Arec))
#####

# Table 2
#####
round(tapply(dataset$DCI,list(dataset$Compet_cod,dataset$V010),mean,na.rm = T),2)

round(tapply(dataset$DCI,list(dataset$Compet_cod,dataset$V010),median,na.rm = T),2)

round(tapply(dataset$DCI,list(dataset$Compet_cod,dataset$V010),sd,na.rm = T),2)

#####

# Figure 1
#####
fig1 <- as.data.frame(round(100*prop.table(table(dataset$V56)),1))

colnames(fig1) <- c("name","Percent")

fig1 %>%
  mutate(name = reorder(name, desc(Percent))) %>%
  ggplot( aes(x=name, y=Percent)) +
  coord_flip()+
  geom_bar(stat="identity") +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(face = "bold", color = "black", size = 16),
    axis.text.y = element_text(face = "bold", color = "black", size = 14))

#####

# Figure 2
#####
fig2 <- rbind(as.data.frame(round(100*prop.table(table(dataset$V57A)),1)),
   as.data.frame(round(100*prop.table(table(dataset$V57B)),1)),
   as.data.frame(round(100*prop.table(table(dataset$V57C)),1)))

fig2$Place <- c(rep("1st place",5),rep("2nd place",5),rep("3rd place",5))

colnames(fig2) <- c("Var","Percent",  "Place")

fig2$Var <- rep(c("At home",
               "City Council",
               "In neighborhood etc",
               "Other than City Council",
               "Other work location"),3)

fig2$Var<-as.factor(fig2$Var)
fig2$Var<-factor(fig2$Var, levels = c("At home", "Other than City Council", "City Council","Other work location","In neighborhood etc"))

fig2$Percent2 <- paste(format(fig2$Percent,digit=2),"%")

#Creating the graph
f2 <- ggplot(data=fig2, aes(x=Var, y=Percent, fill=Place)) +
  geom_bar(stat="identity", position="dodge")

f2 + 
  #coord_flip()+
  #ylab("Percent (%)") + 
  scale_fill_manual(values=c('grey40','gray55','gray70'))+
  geom_text(aes(label=Percent2), vjust=1.6, color="white",position = position_dodge(0.9), size=4.5)+
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(color = "black", size = 16),
    axis.text.y = element_text(face = "bold", color = "black", size = 14),
    legend.position="bottom")

#####

# Figure 3
#####
fig3 <- rbind(as.data.frame(round(100*prop.table(table(dataset$V58A)),1)),
              as.data.frame(round(100*prop.table(table(dataset$V58B)),1)))

# The category "Other" was removed from the graph because it appears just in the first place and is very small
fig3<-fig3[-3,]
###

fig3$Place <- c(rep("1st place",5),rep("2nd place",5))

colnames(fig3) <- c("Var","Percent",  "Place")

fig3$Var <- rep(c("Complaints",
                  "Information",
                  "Requests",
                  "Suggestions",
                  "Report"),2)

fig3$Percent2 <- paste(format(fig3$Percent,digit=2),"%")

#Creating the graph
f3 <- ggplot(data=fig3, aes(x=reorder(Var,-Percent), y=Percent, fill=Place)) +
  geom_bar(stat="identity", position="dodge")

f3 + 
  scale_fill_manual(values=c('gray55','gray70'))+
  geom_text(aes(label=Percent2), vjust=1.6, color="white",position = position_dodge(0.9), size=4.5)+
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(color = "black", size = 16),
    axis.text.y = element_text(face = "bold", color = "black", size = 14),
    legend.position="bottom")

#####

#Figure 4
#####
fig4 <- rbind(as.data.frame(round(100*prop.table(table(dataset$V60A)),1)),
              as.data.frame(round(100*prop.table(table(dataset$V60B)),1)))

# The category "Other" was removed from the graph because it appears just in the first place and is very small
fig4<-fig4[-9,]
###

fig4$Place <- c(rep("1st place",11),rep("2nd place",11))

colnames(fig4) <- c("Var","Percent",  "Place")

fig4$Var <- rep(c("Assistance in doc. acquisition",
                  "Assistance with public services",
                  "Employment requests",
                  "Improvements in the neighborhood",
                  "Legal counsel",
                  "Material inducements",
                  "Medical assistance",
                  "Money",
                  "School admissions",
                  "Support causes",
                  "Transportation"),2)

fig4$Percent2 <- paste(format(fig4$Percent,digit=2),"%")

fig4$Var<-as.factor(fig4$Var)
fig4$Var<-factor(fig4$Var, levels = c("Medical assistance", "Employment requests","Improvements in the neighborhood","Support causes","Material inducements","Assistance with public services","Money","Legal counsel","Assistance in doc. acquisition","Transportation","School admissions"))

#Creating the graph
f4 <- ggplot(data=fig4, aes(x=Var, y=Percent, fill=Place)) +
  geom_bar(stat="identity", position="dodge")

f4 +
  scale_fill_manual(values=c('gray55','gray70'))+
  geom_text(aes(label=Percent2), vjust=1.6, color="white",fontface = c("bold"), position = position_dodge(0.9), size=3)+
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(color = "black", size = 8,angle = 60,vjust = .95,hjust = 1),
    axis.text.y = element_text(face = "bold", color = "black", size = 12),
    legend.position="bottom")
#####

# Figure 5
#####
fig5 <- rbind(as.data.frame(round(100*prop.table(table(dataset$V61A)),1)),
              as.data.frame(round(100*prop.table(table(dataset$V61B)),1)))

# The category "Other" was removed from the graph because it appears just in the first place and is very small
fig5<-fig5[-2,]
###

fig5$Place <- c(rep("1st place",6),rep("2nd place",6))

colnames(fig5) <- c("Var","Percent",  "Place")

fig5$Var <- rep(c("Help from family and friends",
                  "Refer to agency",
                  "Seek help personnel",
                  "Submit a draft bill",
                  "Submit an petition",
                  "Use of own resources"),2)

fig5$Percent2 <- paste(format(fig5$Percent,digit=2),"%")

fig5$Var<-as.factor(fig5$Var)
fig5$Var<-factor(fig5$Var, levels = c("Submit an petition", "Seek help personnel","Refer to agency","Use of own resources","Submit a draft bill","Help from family and friends"))

#Creating the graph
f5 <- ggplot(data=fig5, aes(x=Var, y=Percent, fill=Place)) +
  geom_bar(stat="identity", position="dodge")

f5 +
  scale_fill_manual(values=c('gray55','gray70'))+
  geom_text(aes(label=Percent2), vjust=1.6, color="white",fontface = c("bold"), position = position_dodge(0.9), size=4.5)+
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(color = "black", size = 12),
    axis.text.y = element_text(face = "bold", color = "black", size = 12),
    legend.position="bottom")
#####

# Figure 6
#####
ggplot(subset(dataset, !is.na(DCI)), aes(x=DCI)) +
  geom_bar(stat="count", fill="blue")+
  theme_minimal() + 
  xlab("DCI") +
  ylab("Frequency") +
  geom_text(aes(y = (..count..),label =  ifelse((..count..)==0,"",scales::percent((..count..)/sum(..count..)))), stat="count",colour="white",vjust = +1.5)

#####

# Figure 7
#####
# Generate the new factor variable "CompetPov" by the two other factor variables "Compet_cod" and "V010"
dataset$CompetPov <- with(dataset, interaction(Compet_cod,V010),drop=TRUE )

dataset$CompetPov <-revalue(dataset$CompetPov, c("No.No"="C_NP", "No.Yes"="C_P","Yes.No"="NC_NP","Yes.Yes"="NC_P"))

dataset$CompetPov <- factor(dataset$CompetPov, levels = c("C_NP", "C_P", "NC_NP","NC_P"))

bwplot(DCI ~ CompetPov, 
       data=dataset,
       aspect = 2, 
       ylab="DCI", 
       xlab="Competition x Poverty")

#####

#Figure 8
#####
f8 <- ggplot(dataset,aes(x=V51, y=DCI)) +       geom_boxplot()+ 
      scale_x_discrete(limits=c("No", "Yes"))+
  labs(x="Support the mayor")
#####

#Figure 9
#####
f9 <- ggplot(dataset,aes(x=V37, y=DCI)) +       geom_boxplot()+ 
  scale_x_discrete(limits=c("No", "Yes"))+
  labs(x="Self-declared right-wing")
#####

#Figure 10
#####
dataset$V009<-as.factor(dataset$V009)
dataset$V009<-factor(dataset$V009, levels = c("Small","Medium","Large"))

f10 <- ggplot(dataset,aes(x=V009, y=DCI)) +       geom_boxplot()+
  labs(x="Self-declared right-wing")
#####

#Figure 11
#####
f11 <- ggplot(dataset,aes(x=V010, y=DCI)) +       geom_boxplot()+
  labs(x="Poverty")+
  scale_x_discrete(labels=c("Non-poor", "Poor"))
#####

#Figure 12
#####
f12 <- ggplot(dataset,aes(x=Compet_cod, y=DCI)) + 
  geom_boxplot()+
  labs(x="Electoral competition")+
  scale_x_discrete(labels=c("High", "Low"))
#####

#### Multilevel Ordinal Model

###Table 3: Likelihood ratio tests of cumulative link models: with and without random effects

dataset$DCI_ord <-as.ordered(dataset$DCI)

# Standard cumulative logit model - no covariates
m00 <- clm(DCI_ord ~ 1 ,data = dataset, link = "logit")

# Random intercept cumulative logit model - no covariates
m0 <- clmm(DCI_ord ~ 1 + (1|V004), data=dataset ,nAGQ=8, Hess = TRUE, link = "logit")

# LRT comparing model m0 (random effects) with model m00 (no random effects)
anova(m0,m00)

### Table 4. Results for three multilevel ordinal models (proportional odds) 

# Complete cases selected. 19 cases were eliminated. 
dataset_c <- complete.cases(dataset)
dataset_c <- dataset[dataset_c,]

## Model 1
m0c <- clmm(DCI_ord ~ 1 + (1|V004), data=dataset_c ,nAGQ=8, Hess = TRUE, link = "logit")

summary(m0c)

# Odds ratios
round(exp(m0c$coefficients),3)

## Model 2
m2c <- clmm(DCI_ord ~ V51 + V37 + (1|V004), data=dataset_c ,nAGQ=8, Hess = TRUE, link = "logit")

summary(m2c)

# Odds ratios
round(exp(m2c$coefficients),2)

## Model 3
m3c <- clmm(DCI_ord ~ V51 + V37 + CompetPov + V009 + (1|V004), data=dataset_c ,nAGQ=8, Hess = TRUE, link = "logit")

summary(m3c)

# Odds ratios
round(exp(m3c$coefficients),2)

